tab.path <- "results/table/"
fig.path <- "results/graph/"
fig.name <- paste0(fig.path, "app_fig_e02.pdf")

yvars <- c("par.onl.sh")

Re$Fyear <- factor(Re$year, levels = c(1913, 1911, 1912, 1914))
fcontrols <- paste0("Fyear*", c(controls))

vars <- c( "Fyear*factor(treat_10)", "Fyear", "factor(treat_10)")

m1  <- felm(as.formula(RegFor( y = yvars[1] , x = c(vars, fcontrols) ,
      FE = "county" , IV="0", clust = "klust" )),
           data= Re)
m2  <- felm(as.formula(RegFor( y = yvars[1] , x = c(vars, fcontrols) ,
      FE = "county" , IV="0", clust = "klust" )),
           data= subset(Re, pop <= 15000))
m3  <- felm(as.formula(RegFor( y = yvars[1] , x = c(vars, fcontrols) ,
      FE = "county" , IV="0", clust = "klust" )),
           data= subset(Re,  (dist.roads <= 2 & pop <= 15000)))


kall <- paste0("Fyear", c(1911,1912,1914), ":factor(treat_10)TRUE")

coefs12 <- c(coefficients(m1)[kall[1]], coefficients(m2)[kall[1]], coefficients(m3)[kall[1]])
coefs13 <- c(coefficients(m1)[kall[2]], coefficients(m2)[kall[2]], coefficients(m3)[kall[2]])
coefs14 <- c(coefficients(m1)[kall[3]], coefficients(m2)[kall[3]], coefficients(m3)[kall[3]])

se12 <- c(m1$cse[kall[1]], m1$cse[kall[1]], m1$cse[kall[1]])
se13 <- c(m1$cse[kall[2]], m2$cse[kall[2]], m3$cse[kall[2]])
se14 <- c(m1$cse[kall[3]], m2$cse[kall[3]], m3$cse[kall[3]])


df <-  data.frame(c(coefs12,coefs13,coefs14), c(se12,se13,se14))
names(df) <- c("coef", "se")
df$year <- c(rep(1911,3), rep(1912,3), rep(1914,3) )
df$Sample <- rep(c("Full Sample", "Pop < 10k", "Roads < 2 km and Pop < 15k"),3)

df$lo10 <- df$coef - 1.65*df$se
df$lo5 <- df$coef - 1.96*df$se
df$up10 <- df$coef + 1.65*df$se
df$up5 <- df$coef + 1.96*df$se

df <- subset(df, Sample != "Pop < 10k")
df <- rbind(df, c( NA, NA, 1913, NA, NA, NA, NA))
kols <- ghibli_palette(name = "PonyoMedium")[c(2,4)]



p.t <- ggplot(data = df, aes(x = factor(year))) +
geom_pointrange(aes( y= coef, ymin=lo5, ymax = up5, col = Sample, shape = Sample) , 
      size = 2, fatten = 1,  alpha = 0.9, position = position_dodge(width = 0.4)) +
geom_linerange(aes( ymin=lo10, ymax = up10, col = Sample) , 
      size = 3,  alpha = 0.35, position = position_dodge(width = 0.4)) +
scale_color_manual(values = kols, na.translate=FALSE) +
scale_shape_manual(values = c(15, 18, 19),na.translate = FALSE) +
geom_vline(xintercept=3,linetype="dashed", size=0.2) +
geom_hline(yintercept=0,linetype="dotted", size=0.2) +
theme_light(base_size = 16) + labs( x = "Year" , y = "Effect Year x March") +
theme(legend.position = "bottom",
        legend.direction = "vertical") +
guides(shape = guide_legend(override.aes = list(size = 0.5)))

ggsave(fig.name, plot = p.t, height = 5.6, width = 8 )
rm(list=
    setdiff(ls()[sapply(ls(), function(x) any(class(get(x)) == 'data.frame'))],
    c("Re","Ind")))